library(Seurat)
library(SeuratData)
library(ggplot2)
library(patchwork)
library(dplyr)

Loading data

InstallData("stxBrain")

brain <- LoadData("stxBrain", type = "anterior1")

Data preprocessing

plot1 <- VlnPlot(brain, features = "nCount_Spatial", pt.size = 0.1) + NoLegend()
plot2 <- SpatialFeaturePlot(brain, features = "nCount_Spatial") + theme(legend.position = "right")
wrap_plots(plot1, plot2)

brain <- SCTransform(brain, assay = "Spatial", verbose = FALSE)

Gene expression visualization

HPCA is a marker for hippocampus tissue, and TTR is a marker for the choroid plexus.

SpatialFeaturePlot(brain, features = c("Hpca", "Ttr"))

# some parameters in this function: pt.size.factor (default=1.6) = the size of the spots displayed 
# alpha (default c(1,1)) = min and max transparency. as expression levels increase, transparency increases. 
p1 <- SpatialFeaturePlot(brain, features = "Ttr", pt.size.factor = 1)
p2 <- SpatialFeaturePlot(brain, features = "Ttr", alpha = c(0.1, 1))
p1 + p2

Dimensionality reduction, clustering and visualization

DimPlot() plots the UMAP and SpatialDimPlot() overlays the clusters on the image of the tissue.

brain <- RunPCA(brain, assay = "SCT", verbose = FALSE)
brain <- FindNeighbors(brain, reduction = "pca", dims = 1:30)
brain <- FindClusters(brain, verbose = FALSE)
brain <- RunUMAP(brain, reduction = "pca", dims = 1:30)

p1 <- DimPlot(brain, reduction = "umap", label = TRUE)
p2 <- SpatialDimPlot(brain, label = TRUE, label.size = 3)
p1 + p2

Can also use the cells.highlight parameter to visualize the location of individual clusters on the image in the SpatialDimPlot function:

SpatialDimPlot(brain, cells.highlight = CellsByIdentities(object = brain, idents = c(2, 1, 4, 3,
    5, 6)), facet.highlight = TRUE, ncol = 3)

Can also makes these plots interactive. More information in the vignette.

Identification of spatially variable features

There are 2 workflows to identify molecular features corresponding to spatial locations on tissues. The first works well in cases where anatomical annotations are available for the tissues. This works well for this dataset since the different regions of the brain are prior knowledge here, and the clustering aligns with the different anatomical regions. The results of this method are shown below.

de_markers <- FindMarkers(brain, ident.1 = 5, ident.2 = 6)
SpatialFeaturePlot(object = brain, features = rownames(de_markers)[1:3], alpha = c(0.1, 1), ncol = 3)

The second method looks for molecular features that show spatial patterns on the tissue when anatomical annotations are not available. This is accomplished through the FindSpatiallyVariableFeatures function.

In the vignette, they subset the top 6 spatially variable features and plot their distribution on the tissue. However, this line of code gave an error. To solve this, I found a solution that re-writes the function, named SpatiallyVariableFeatures_workaround(). This worked as expected. The results are shown below.

brain <- FindSpatiallyVariableFeatures(brain, assay = "SCT", features = VariableFeatures(brain)[1:1000],
    selection.method = "moransi")

# top.features <- head(SpatiallyVariableFeatures(brain, method = "moransi"), 6)

SpatiallyVariableFeatures_workaround <- function(object, assay="SCT", selection.method = "moransi") {
  #' This is work around function to replace SeuratObject::SpatiallyVariableFeatures function.
  #' return ranked list of Spatially Variable Features
  
  # Check if object is a Seurat object
  if (!inherits(object, "Seurat")) {
    stop("object must be a Seurat object")
  }

  # Check if assay is a valid assay
  if (!assay %in% names(object@assays)) {
    stop("assay must be a valid assay")
  }
  
  # Extract meta.features from the specified object and assay
  data <- object@assays[[assay]]@meta.features
  
  # Select columns starting with the provided col_prefix
  moransi_cols <- grep(paste0("^", selection.method), colnames(data), value = TRUE)

  # Filter rows where "moransi.spatially.variable" is TRUE
  filtered_data <- data[data[[paste0(selection.method, ".spatially.variable")]], moransi_cols]

  # Sort filtered data by "moransi.spatially.variable.rank" column in ascending order
  sorted_data <- filtered_data[order(filtered_data[[paste0(selection.method, ".spatially.variable.rank")]]), ]

  # Return row names of the sorted data frame
  rownames(sorted_data)
}

top.features <- head(SpatiallyVariableFeatures_workaround(brain, selection.method = "moransi"), 6)

SpatialFeaturePlot(brain, features = top.features, ncol = 3, alpha = c(0.1, 1))

Subset out anatomical regions

Next, they demonstrate how to isolate spots/cells from a specific anatomical region of the tissue. In this case, they subset only those spots/cells that belong to clusters that represent the cortex region.

I was able to create the initial cortex Seurat object as a subset of the clusters. However, I was unable to further refine the cell/spot selection to just the cortical region using the code that was provided. It looks like variable names are incorrect. Specifically, when I looked at the images layer in the cortex object, I was unable to find variables “anterior1_imagerow” and “anterior1_imagecol”.

cortex <- subset(brain, idents = c(1, 2, 3, 4, 6, 7))
# now remove additional cells, use SpatialDimPlots to visualize what to remove
# SpatialDimPlot(cortex,cells.highlight = WhichCells(cortex, expression = image_imagerow > 400 | image_imagecol < 150))

#cortex <- subset(cortex, anterior1_imagerow > 400 | anterior1_imagecol < 150, invert = TRUE)
#cortex <- subset(cortex, anterior1_imagerow > 275 & anterior1_imagecol > 370, invert = TRUE)
#cortex <- subset(cortex, anterior1_imagerow > 250 & anterior1_imagecol > 440, invert = TRUE)
p1 <- SpatialDimPlot(cortex, crop = TRUE, label = TRUE)
p2 <- SpatialDimPlot(cortex, crop = FALSE, label = TRUE, pt.size.factor = 1, label.size = 3)
p1 + p2